Efficient Computation of Closed-loop Frequency 

Response for Large Order Flexible Systems 
Peiman G. Maghami 
NASA Langley Research Center 
Daniel P. Giesy^ 

Lockheed Martin Engineering & Sciences 

Abstract 

An efficient and robust computational scheme is given for the calculation of the 
frequency response function of a large order, flexible system implemented with a 
linear, time invariant control system. Advantage is taken of the highly structured 
sparsity of the system matrix of the plant based on a model of the structure using 
normal mode coordinates. The computational time per frequency point of the new 
computational scheme is a linear function of system size, a significant improvement 
over traditional, full-matrix techniques whose computational times per frequency 
point range from quadratic to cubic functions of system size. This permits the 
practical frequency domain analysis of systems of much larger order than by 
traditional, full-matrix techniques. Formulations are given for both open and closed 
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loop systems. Numerical examples are presented showing the advantages of the 
present formulation over traditional approaches, both in speed and in accuracy. 
Using a model with 703 structural modes, a speed-up of almost two orders of 
magnitude was observed while accuracy improved by up to 5 decimal places. 

Introduction 

Control of flexible systems has received significant attention in the literature. 
To date, numerous techniques, algorithms and procedures have been developed for 
design of controllers for such systems ranging from spacecraft and satellites to 
aircraft, ships, machines, etc. These flexible systems which are generally infinite- 
dimensional are typically modeled using a finite number of generalized coordinates 
or modes. Control of flexible systems may become difficult depending on the 
number, location, relative proximity, and inherent damping of these modes. The 
response of the system to a given disturbance/excitation generally depends on 
modal properties (amplitude, frequency, and damping) and the amplitude and phase 
content of the disturbance/excitation. In general, two techniques, time domain 
analysis and frequency domain analysis, have been developed and extensively 
used to analyze and characterize the input/output behavior of linear time-invariant 
systems including flexible systems. In frequency domain analysis, frequency 
response functions (defined as transfer function matrices from inputs to outputs of 
the system) have typically been used (usually in the form of magnitude and phase 
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or Bode plots) in the analysis of linear systems as well as in designing controllers 
for such systems. In general, frequency response functions of the open-loop system 
are used to evaluate the performance of the open-loop system, and to identify and 
quantify needed performance and/or stability improvements at various frequency 
bands. The closed-loop frequency response functions are typically needed to insure 
that desired performance and stability have been achieved by the control system. 
Moreover, frequency-domain specifications such as peak magnitude, bandwidth, 
roll-off rate, and etc. are often used in characterizing the desired behavior of the 
system in the frequency domain (this is known as loop shaping). 

In general, the order of the flexible system (as defined by the number of modes 
retained in the model) for which open-loop and/or closed-loop analysis is performed 
depends on the application considered. For example, if the closed-loop response 
of a spacecraft with a low-bandwidth attitude control system is of interest, then a 
small set of modes would be sufficient to capture the low frequency closed-loop 
behavior of the system. On the other hand, if the response of the flexible system 
is desired over a large frequency range or if the control system considered has a 
high bandwidth, then a large set of modes (in the hundreds or thousands) may be 
necessary to capture the true response of the system. 

However, the current techniques for obtaining frequency response functions, 
although able to deal with small or medium size systems, have problems in handling 
large order systems. A straightforward calculation of the frequency response 
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function matrix at a single frequency point which is based on the definition of 
the transfer function has a computational cost which is a cubic function of the 
system size. If this calculation must be repeated for many frequency points, Laub 
([1] and [2]) presents a technique which has a better average cost. This technique 
performs an initial orthogonal transformation of the system which reduces the 
system response matrix to Hessenberg form. This initial transformation has a 
computational cost which is a cubic function of the system size. This technique 
can then calculate the frequency response function matrix at each frequency point 
at a cost which is a quadratic function of the system size. However, for very large 
systems (many hundreds of modes or more), even this is too slow, and a better 
method is needed. 

To this end, this paper describes a novel and efficient technique for the 
computation of closed-loop frequency response functions of large order flexible 
systems. The proposed technique is computationally robust and accurate. It 
takes advantage of the sparsity of the flexible systems in normal mode coordinates 
and reduces the computational cost from a quadratic function of the order of the 
system to a linear function. Formulations are given for both open and closed loop 
systems. Numerical examples are presented showing the advantages of the present 
formulation over traditional approaches, both in speed and in accuracy. 
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Mathematical Formulation 


Second-Order Modal Equations 

The dynamics of a typical linear, time-invariant flexible system may be written 
in a second-order form as 

Mx + D.r + Kx = Hu + Hrfic 
y = C p x + C r x 
}l pr — C'P r r I C'P r T X- C'P r T 

(-J V / p kIj | E/ 1 | v_x if 


where AT, D, and A' are the n x n mass, damping and stiffness matrices, respec- 
tively; x is the n x 1 position/attitude vector; // is the in x 1 control input vector; w 
is the r x 1 disturbance vector; H is the n x m control input influence matrix; and 
H ( i is the n x r disturbance influence matrix. The vectors // and are, respec- 
tively, the q x 1 measurements output vector and the / x 1 vector of performance 
outputs; Cp and C r are q x n measurement output influence matrices; and Cjf, 
Cf‘, and cr are / x n performance output influence matrices. 

If the second-order system is transformed into normal mode coordinates, and p 
of the normal modes are retained to capture the relevant dynamics of the structure, 
then the system equations may be written in a modal form as 

Mi) | Dq + Kq = Hu + H ( /ir 
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y = C p q + C r q 


y pr = Cfq + C pr q + C pr q 


where M, D, and K are the p x p modal mass, damping, and stiffness matrices, 
respectively; q is the px 1 vector of modal coordinates; and H and Hy are the pxm 
control input and the p x r disturbance influence matrices in modal coordinates, 
respectively. The matrices C p and C r are <j x /> measurement output influence 
matrices in modal coordinates; and C pi , C pi , and C pi are / x /> performance 
output influence matrices in modal coordinates. 

It is assumed that the mode shapes are normalized with respect to the mass 
matrix, and modal damping is assumed. This means that M = I p , D = 
diag{2(ia;i, 2 ( 2 ^ 2 , ...,2^^}, and K = diag{u;f, uj\, ...,u) p \ where I p is 
the identity matrix of order p and up and are the open-loop frequencies and 
damping ratios. 

The control input and disturbance influence matrices are given by: 

H = § t H 

Hi = $ T H d 

The measurement and performance output influence matrices are given by: 

C p = C p $ ; C r = C r fl> 
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(f = Cf<5> ; Cf' = ; Cf = 


The columns of matrix <f> are the p retained mode shapes: 


d> = [0i 02 • • • 0p . 


The second-order modal equations may be rewritten in a first-order form as 

x s = • fs.r.s + + B C [W 


y = C'xs (l) 

yP r _ (jF> Xg _|_ C ^ 1 X s . 

The vector ./\ s is the plant state vector whose components are 

'qi ' 

01 
Q2 

j & 

X s = < . 


Qp 
Qp 

and the vectors // and ij pr are the same plant measurement and performance outputs, 
respectively. The matrix A s is the plant state matrix and has the form 


r 4 0 


A, 


0 


0 A 2 ■■■ 0 


where 


A[ 


0 0 


0 


A* 


(2) 


1 


-LO: 


( 3 ) 


7 



The matrix B s is the control input influence matrix, formed by setting its odd- 
numbered rows to zeros and using the rows of H for its even-numbered rows: 

0 0 0 

#11 #12 #lm 

0 0 0 

#21 #22 # 2 ??? 

0 0 0 

Hp\ Hp2 Hpm 

in which Hjj, for example, represents the ( i. j ) element of matrix H. The matrix 
B ( i is formed from Hq in the same manner. 

The measurement output influence matrix, C is defined by setting the odd- 
numbered columns of C to the columns of C p and the even numbered columns 
of C to the columns of C r \ 

'<5,(1, 1) <5 r (l,l) <5,(1, 2) <5 r (l,2) <5,(1, p) C r ( l,p) 

<5,( 2,i) C'r(2,l) <5,( 2,2) C'r(2,2) <5,(2, p) C r (2,p) 

C= ; ; 

Cp{ c L 1) Cr(q, 1) Cp{q, 2) C' r (q, 2) Cp(q,p ) C r (q,p ) 

where C p ( /. /) and C r (i,j) denote the (i,j) element of matrix C p and C r , 
respectively. C l " is defined from Cp 1 and C 1 *' in the same fashion. CfJ' is defined 
by setting the odd numbered columns of to zeros and the even numbered 
columns of C f ? to the columns of CJfJ . 




By substituting the first equation of (1) into the third, the acceleration term can 



be replaced by feedthrough: 

d' s = A s x s + B s u + Byw 

y = Cx s (5) 

yP r = C pr Xs + ]jpr u + £>£> 

The performance output influence matrix is given by 

C* r = Cf + cf .4, 


while the performance feedthrough matrices are 

D 1 " = CfB a ; D 1 " = CfB d . 

Notice that if there is no performance acceleration output (Cf ? = 0), then 
Cl r = 0 and C* r = 0, so both feedthrough matrices, Du and D 1 !,' , are zero. 

Control System Equations 


In this paper, it is assumed that the structure is controlled by a linear time- 
invariant control system. The model of a linear time-invariant control system for 

a typical flexible structure may be written as 

x c = A c x c + B c y 


u = —C c x c 

where x c denotes the /.■ x 1 vector of control system states; A c , B c , and C c represent 
the k x k control system state matrix, the k x q input influence matrix, and the 
m x k output influence matrix, respectively; and y is the measurement output vector 
which was defined in the previous section. 
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Frequency Domain Equations 


For the open-loop plant, y and u of equation (5) are nonexistent, so the equations 
reduce to 


is = A s x s + Byw 
y pr = C pr x s + Dfiw. 


The open-loop transfer function from the disturbances, w, to the performance 
output, y pr , is defined for all complex s not in the spectrum of A s and is given by 

T(s) = C' pr {sl - A s )~ 1 B d + D%. (7) 


The closed loop system is more complicated. Using equations (5) and (6), the 
closed-loop dynamics of the controlled structure may be written as 

x = .4a: + Byw 


( 8 ) 


y 


pr = C pr x + DP w r 


w 


where x represents the closed-loop state vector defined as 


x = 



(9) 


.4,5, and CJ pr are, respectively, the closed-loop state matrix, disturbance influence 
matrix, and output influence matrix. These matrices are given by: 


.4 


A s —B s Cc 
B c C A c 


; By 


By 

0 


C pr = [C pr —D pr C c ] 


( 10 ) 


to 



The closed-loop transfer function from the disturbances, w, to the performance 
output, y pr , is defined for all complex s not in the spectrum of A and is given by 

f(s)=C‘’ r (.sI-Ay l B d + D>". (11) 


Open-Loop Calculation 


The algorithm presented here for calculation of the frequency response function 
of an open-loop structural system seems to be a part of engineering folklore. It is 
presented here for completeness and because it is a building block for the closed- 
loop algorithm to follow. 

Assume s is not in the spectrum of A s . From equations (2) and (3), it follows 
that (si — As) -1 is block diagonal with the i-th block being 


(A - 4) 



+ C, ; — •; -s 


Aujf 


s A 2 Cicj-i 




i = 1,2,. . . ,p; 


where I-j denotes the 2 by 2 identity matrix. Furthermore (see the discussion 
following equation (4)) the odd-numbered rows of B c i are zero. If the row i of 

( ,j\ -j 

B ( i is denoted by b (l , if Q(s) is used to represent (si — A s ) B ( i, and if (5(s)is 
partitioned as 


Q(s) 


Q i(s) 
Q'2(s) 

_Qp ( 5 ) _ 


( 12 ) 


li 



where each partition matrix Qi{s) is a 2 by r matrix, then Q is calculated using 
the formula 


Qi{s) — ^1/ ^5 2 + 2QuJiS + 




; * = 1,2,.. Mp. 


(13) 


The open-loop transfer function calculation is completed in a straightforward 


manner: 


T{s) = Cr r Q(s)+D% 


If desired, the gain and phase angle data (Bode plot data) may then be computed 
directly from the frequency response function matrix. If there are no acceleration 
performance measurements, then the feed-forward term is zero and the software 
may bypass the step where D 1 " is added in. 

The standard, full matrix way to calculate Q{s) would involve first performing 
an LU decomposition of si — A s followed by a backward and then forward solve 
of the triangular systems of equations using the columns of B ( i as right-hand sides. 
The FLOP (FLoating point Operations) count for this is 0(p 3 ) + 0(p 2 r), so T(s) 
is computed in 0(p 3 ) + 0(p 2 r) -\-0{plr) FLOPS. Thus, in the typical case where 
system size is much larger than the number of disturbances, the calculation time 
per frequency point is a cubic function of system size. 

If this calculation must be repeated for many values of s (a typical scenario), 
the technique of [1] and [2] has a better average FLOP count. An initial 0(p 3 ) 
orthogonal transformation must be done once; then for each s, Q(s ) is calculated 
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in 0{p 2 r) FLOPs, so T(s) is computed in 0{p 2 r ) + O(plr) FLOPS. Thus, if 
the number of frequencies for which this calculation must be repeated is on the 
order of the system size, or larger, the calculation time per frequency point is a 
quadratic function of system size. 

When the calculation is done as in equations (12) and (13), the flop count is 
O(pr), so T(s) is computed in O(plr) FLOPS. Thus, the calculation time per 
frequency point is a linear function of system size. This represents a substantial 
savings, particularly when a large number of modes is necessary to capture the 
dynamics of the system. 

Closed-loop Calculation 

The closed-loop dynamics of the controlled system are given in equations (8), 
(9), and (10). 

Observing the closed-loop state matrix A, it is obvious the block diagonal 
form of the open-loop plant has been destroyed by the coupling generated by the 
feedback connection of plant and the control system. However, the initial sparsity 
of the open-loop state matrix A s is still intact. Now, the sparsity of the open-loop 
state matrix is exploited to develop an efficient method for the computation of 
closed-loop frequency response function matrix of the controlled flexible structure. 
If sparsity is not exploited and many structural modes are modeled, it follows 
from equation (11) that a large computational effort would be required to calculate 
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the closed-loop frequency response function matrix, since this would involve the 
computation of the matrix term (si — Aj B ( i with s = juj for all desired 
frequency values uj. 

In the following it is assumed that s is not in the spectrum of A (necessary for 
the transfer function even to be defined) and it is further assumed that s is not in 
the spectrum of A s . This further assumption is needed to enable some algebraic 
manipulation, and should not adversely affect the applicability of the following 
results. On the one hand, since A s is the plant state matrix for a linear model of 
a flexible structure, its eigenvalues occur either at 0 (corresponding to rigid body 
modes) or in the left half plane (corresponding to damped flexible modes). On the 
other hand, it is anticipated that these results will be used to compute T(s) for 
s = juj with to > 0. Thus, excluding the eigenvalues of A s from the domain of 
applicability of these results does not impact the anticipated usage. 

The matrix term fsT — A\ in equation (11) may be written as 


AA' 

l '~H 

1 

aAA 

II 

slg Ag 

BgC' c " 

= 

'En(s) 

E\2 

-B c C 

sl c Ac 


E-21 

E-22{s) _ 


where I s , and I c are identity matrices of orders equal to the size of plant state 
vector and controller state vector, respectively. Introduce the notation: 


'A'n(s) 

Ai 2 (s)' 

= ( s j — 4^ 1 — 

'En(s) 

El2 

X 2 i( S ) 

A 22 (s)_ 

- V 1 A ) - 

E-2 1 

E22{s) _ 


( 15 ) 
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The assumptions which have been made about s insure that the inverses in (15) 


exist, as does . Rewrite (15) as: 

A'n A r i2 Eh E12 _ Is 0 

A 2 1 A 22 J \_E '21 E '22 \ |_ ° Ic\ { } 

Expanding the lower left block of (16) and solving for A'21 gives A'21 = 
— A 22 1 AP| 1 . Expanding the lower right block of (16), substituting the previous 

expression for A'21, and factoring gives A22 A = I c , where A = E 22 —E 21 EjQ 1 E\ 2 . 
This demonstrates that A is invertible, and justifies the application of the block 
matrix inversion formula given in [3, page 898] to find the inverse of the block 
matrix in Equation (15): 

A = E22 — E21E11 E12 
A'n = E11 + E11E12 A 1 i?2i-E'n 1 
A 12 = -E^EuA- 1 (17) 

A'21 = —A 1 E2iEii 

X 22 = A ” 1 

Using equations (10) and (15) in equation (11), the closed-loop transfer function 
from the disturbances to the performance outputs is reduced to: 

T( S ) = [o ir -Di r c c ] B 0 d +d>;; 

= C» r XnB d - DflaXnB,, + D!”’ 

Using equation (17) this becomes: 

f =C 1>r Eli B c [ + C‘ ,r El 1 1 E u A- 1 E-nEl 1 1 B i 
+ N’l'C'cA - 1 E21EI1 B d + D% 
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Using equation (14), replace £12 and E-)\ in: 


1. the expression for A in equation (17) 

2. the previous equation 

This produces: 

A =E-22 + B c CE u B s C c 

T =C 1>V Eyy Bd - Cl’ r E^B.,{c c A- 1 B c CE^B d '} (18) 

- Dl’ {c c A- 1 B c CE^B d } + D £, r 

In this form, the following computational efficiencies are observed: 

□ Since E\ \ = (si — A s ) and B s shares with B ( i the property of having zeros in 
the odd numbered rows, the terms E^Bs and Cf/ B ( i can be computed using 
the techniques presented for efficient computation of the open-loop transfer 
function (see equation (13)). 

□ The computation of A ~ 1 B C must be done as a full matrix computation, but 
since A is of the same order as the control system, which is usually small 
compared to the order of the analysis model of the plant, it should not be very 
costly to compute. 

□ Common sub-expressions, such as those mentioned in the previous items and 
those enclosed in braces in equation (18) are computed once per frequency, 
saved, and reused. 

□ The expected shapes of the matrices and the exploitation of common sub- 
expressions make it advisable not to precompute some of the matrix products 
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in equation (18) which are independent of the frequency parameter s (if X. is 
a tall, skinny matrix, and Y is a short, wide matrix, so that W = XV is both 
tall and wide; and if is a vector, then calculating X(Yz) is cheaper than 
calculating Wz). 

□ If there are no acceleration sensors in the performance outputs, so that the 
feed-forward matrices are zero, the software may bypass the second line of the 
T calculation in equation (18). 

Now, using equation (18) to calculate T(s), the frequency response function 
matrix of the closed-loop system is evaluated for various values of s = juj, with 
uj taking on the user-specified frequency values. The closed-loop gain and phase 
plots (Bode plots) may then be computed directly from the frequency response 
function matrix, if desired. 

The closed loop system matrix A (equation (10)) has order 2 p + k. As in 
the discussion of the open loop calculation, if T(s) is calculated as in equation 
(11) using standard full matrix techniques, the computation takes o{[2pY k)^j + 
o(^(2pY k) 2 r^j + ()( (2/> + k)lr) FLOPs per frequency point. Using the technique 
of [1] and [2], if the number of frequency points for which the calculation is 
to be repeated is on the order of 2 p + k or more, then T(s) can be calculated 
in O (^(2p A k) 2 r^j A ()((2p + k)lr) FLOPs per frequency point. Again, these 
are cubic and quadratic functions, respectively, of the system size. By counting 
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FLOPs resulting from subroutine calls and DO loops in the FORTRAN software 
used to implement the calculation of equation (18), it is determined that T(s) can 
be calculated in 

0(p(lr + rq + qm + rm)) + O^A : 3 + k 2 r + k(rq + qm + rm)j + O(lrm) 

FLOPs per frequency point. Again, this is a linear function of system size. 

In future work, it is expected that the 0{k 3 ) in this last expression, which 
comes from performing a LU decomposition on the matrix A, will be reduced to 
0(k 2 ). This will be based on applying the technique of [1] and [2] to E 22 and 
making use of the observation that the other term in the definition of A is, in the 
expected applications, of low rank. 

Software Implementation 

The evaluation of the open- and closed-loop transfer function has been imple- 
mented using MATLAB function M-file (MATLAB, a product of The Math Works, 
Inc., is “a technical computing environment for high-performance numeric com- 
putation and visualization” [4, page /]), and as FORTRAN 77 code which is then 
accessed through MATLAB using the MEX-file external interface facility. The 
M-files contain straightforward implementations of the calculations presented in 
the preceding two subsections. 

The FORTRAN source code for the MEX-files uses the Basic Linear Algebra 
Subprograms (BLAS, [5], [6], [7], [8], [9], [10], and [11]) to perform vector- 
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vector, vector-matrix, and matrix-matrix operations. In addition, LAPACK ([12]) 
subroutine ZGESV, a complex double precision linear equation solver, is used to 
calculate A ~ 1 B C . 

Numerical Examples 

A number of numerical examples are presented to demonstrate the efficiency 
and accuracy of the algorithm presented in this paper compared to two standard 
full matrix methods of calculating the frequency response function of a closed- 
loop system. 

The EOS-AM-1 Spacecraft Model 

The data comes from a model for the EOS-AM-1 spacecraft used in a jitter 
reduction study ([13]). The structural model contains 703 modes for a potential 
1406 plant states. There are 6 rigid body modes and flexible modes ranging 
from 1.24 to 1564 radians per second. The 6 measurement outputs are the 
spacecraft’s roll, roll rate, pitch, pitch rate, yaw, and yaw rate measurements at 
the spacecraft navigational unit. Actuators consist of x-, y-, and z-axis torquers. 
The control system has 39 states. Up to 10 channels of disturbance input and 
27 channels of performance measurement output were used. Each case was run 
using position measurements at each output, resulting in no feedforward term, and 
using acceleration measurements at each output, resulting in a feedforward term 
being present. 
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All algorithms used in this timing study are intended to be used to calculate 
the frequency response matrix at multiple frequency points, so that the frequency 
response may be plotted (as, e.g., Bode plots). They all have some calculations 
which are done once per entry into the algorithm and other calculations which are 
done once for every frequency point. To take account of this, all cases were run 
over a range of 200 frequency points and most of them were rerun using 2000 
points (the exceptions were the cases which would have required 5+ days of cpu 
time to complete). In all cases, the points were logarithmically distributed between 
frequencies of .01 and 10000 radians per second. 

Software Used in Timing Studies 

In this study, two software realizations of the closed-loop frequency response 
function calculations are compared to two software realizations of previously 
available algorithms. 

The present algorithm is programmed both as a MATLAB function M-file and 
as FORTRAN 77 code which is then accessed through MATLAB using the MEX- 
file external interface facility. These will be called, respectively, the new M-code 
and the new Mex-code. 

One of the programs used for comparison makes use of the algorithm in [1] 
and [2]. The FORTRAN code in [2] is in single precision; Laub’s own double 
precision FORTRAN code is imbedded in the software package FREQ ([14]) and 
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was used here. This test code is purely in FORTRAN 77. This will be called the 
old FORTRAN code. 


Preliminary testing indicated that, in order to get a reasonably well-conditioned 
matrix for the si — A expression in the Laub code, it was necessary to exercise the 
built-in option of balancing the A matrix. The unbalanced matrix was particularly 
ill-conditioned at low frequencies. This can be attributed to the presence of the 
rigid body (0 frequency) modes. In the Laub code, balancing was coupled with 
the extraction of the eigenvalues of A. As Laub wrote this code, the same value 
of the input flag which signaled the code to balance the A matrix also signaled the 
code to extract its eigenvalues. For purposes of timing tests here, the Laub code 
was modified so that the portion which extracts eigenvalues was bypassed. 

The other program used for comparison is the Math Works M-file freqrc.m, 
an undocumented utility routine in the Robust Control Toolbox, [15], which calcu- 
lates (to quote the program preamble comments) “Continuous complex frequency 
response (MIMO)”. This will be called the old M-code. Once again, to achieve 
reasonable accuracy, it was necessary to balance A. This was done using MAT- 
LAB built-in routine balance. 

Timing Comparisons 

The executions times of the four test codes are compared on 12 representative 
problems. Three different plant sizes were used: a small plant with 1 input, 1 
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output, and 61 states (24 structural and 39 controller states); a medium plant with 
3 inputs, 5 outputs, and 221 states (184 structural and 39 controller states); and 
a large plant with 10 inputs, 27 outputs, and 1445 states (1406 structural and 
39 controller states). For each plant size, two plants were used: one with no 
feedforward term (i.e., no acceleration outputs were used as performance outputs) 
and one with feedforward. The frequency response function of each of these 6 
plants was calculated at a short (200 values) vector of frequency values and at a 
long (2000 values) vector of frequency values. 

Particularly on the larger plants, the new algorithm performs dramatically faster 
than the older programs. This should not be taken as an indictment of the older 
technology. The older technology was designed to apply to an arbitrary plant while 
the new takes full advantage of the particular pattern of sparsity which results from 
using the modal model of a flexible structure. On the other hand, when the new 
technology is applicable, it enables analysis of structures of much larger order than 
would be practical or even possible with the older technology. 

Table 1 gives the time in seconds to calculate the frequency response function 
using each of the 4 test routines for each of these 12 cases (except that the old 
M-code does not attempt the two largest cases). 

One conclusion to be drawn from this table is that the timing values returned 
by the system timing software are not totally consistent with each other. The first 
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three software packages in that table all check up front to see if feedforward is 
present. The bulk of the code is executed whether feedforward is present or not. 
If feedforward is present, additional code is executed which should take additional 
time. But in 9 of 18 cases, the table shows the feedforward case taking less time 
than the one without. 

That said, there are still significant trends to be observed in this timing data. The 
new Mex-code is significantly faster than the M-code; in the more important larger 
cases, about 3 times as fast. This justifies the effort of rendering the algorithm in 
FORTRAN and writing the interface necessary to access it through MATLAB. 

Comparing the new Mex-code, which is FORTRAN based, with the old FOR- 
TRAN code shows that for the small system, the old code more than holds its 
own. This is not unexpected, since in the small system, the controller dominates 
the count of states. Thus, there is relatively little of the sparsity from the structural 
part of the plant of which the new Mex-code may take advantage. But even in the 
medium size case, the new code is 4 to 7 times as fast as the old. This is getting 
near the size at which conventional numerical analytic wisdom would place the 
limits of applicability of the old, full matrix based, technique. Since the time for 
the old FORTRAN code is expected to grow quadratically with the number of sys- 
tem states while that of the new technique is expected to grow only linearly, it is 
not surprising that the difference between them in the largest example is so great. 
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In what is called here the old M-code, Math Works actually used M-hles for 
outer loop logic control to drive a built in function, itifr. This function calculates 
the matrix G whose columns are (s(i)I — A)~ 1 b where s is a vector of complex 
numbers (set in the present application to ju, where j 2 = — 1 and cc is a vector 
of frequencies) and b is a column vector (set in this application to a column of the 
B matrix). The on-line help for itifr states that it “implements, in high speed” 
what the user could calculate by looping through the elements of s and building G 
one column at a time. Despite this, it is only competitive on the smallest system, 
and then only against the new M-code which utilizes MATLAB built-in functions 
only at the more primitive level of basic matrix operations. 

All of the algorithms tested do have some “once per entry” calculations in 
addition to the calculations which occur once per frequency value. Thus, the time 
for the 2000 point calculations should never be more than 10 times that for the 
200 point calculations. In Table 1, there are several exceptions to this. This 
reinforces the previous remark that the numbers returned by the computer system 
timing routines are, at best, approximate. However, from looking at the largest 
case, it can be reasonably concluded that the “once per entry” overhead is fairly 
small in both realizations of the new algorithm while being substantial in the old 
FORTRAN code, at least for large systems. This is expected, since for a system of 
order n (all other parameters being held fixed), the “once per entry” overhead in 
the old FORTRAN code includes the initial reduction which takes 0(n 3 ) FLOPs, 
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while the “per frequency point” calculation takes 0(n 2 ) FLOPs. 

Accuracy 

No formal error analysis has been performed on the new algorithm. There is, 
however, numerical evidence to support the thesis that the new algorithm is more 
accurate than the older techniques, particularly when applied to larger systems. 

Outputs from the four algorithm realizations were compared. For each fre- 
quency value, individual entries in the frequency response matrices computed by 
the four codes were compared using a symmetric relative error: The discrepancy 
between complex numbers and w (not both 0) was measured by 
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This error measure ranges from a minimum of 0 to a maximum of 2. A value 
of 6(z,w) near HP" indicates that and w agree to about n decimal places 
while S(z,w) > .1 indicates anything from rough approximation (near .1) to no 
correlation (bigger than, say, 1). For each fixed frequency, the worst discrepancy 
over all possible input-output pairs was observed. 

The size of the discrepancy between the frequency response function matrices 
computed by these codes was observed to depend not only on which two of the 
codes were being compared but also on the size of the system, the frequency, and 
whether or not feedforward was present. It would take too much space to present 
details of these comparisons. However, some general statements can be made. 
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For each of the test problems (corresponding to one row of Table 1), the results 
produced by the four codes were compared two by two. The overall best agreement 
between any pair of calculations came from comparing the outputs of the new 
Mex-code and the new M-code. At worst, these agree to about 7 decimal places. 
This generally improves with reduction of system size or increase in frequency so 
that best agreement is within machine accuracy. No other pairing, either between 
one of the new codes and one of the old or between the two old, ever showed 
noticeably better agreement, and in general the agreement was much worse. It 
frequently occurred that the results from comparing the two new codes showed 
that the agreement of their computations was better than that of any other pairing 
by at least 2 decimal places. In the largest system, this advantage could increase 
to 5 decimal places, particularly for small frequencies or when no feedforward 
was present. 

Thus, two dissimilar implementations of the new algorithm produce results in 
good agreement. When a parallel process is applied to the older algorithm, the two 
dissimilar implementations produce results which are not in such good agreement, 
either with each other or with those of the new algorithm. 

To provide further evidence that the results of the new algorithm are the more 
accurate, the old FORTRAN code was translated to quadruple precision from its 
native double precision and was run (at a time penalty of about 32 x ) on the medium 
sized problem using no feedforward and 200 frequency points. The output from 
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this showed the same degree of agreement with the output from the new codes as 
they showed with each other. 

These results combine to indicate strongly that the new algorithm provides 
more accurate results than those previously available. There are theoretical grounds 
for expecting this. The old way requires the solution of linear systems with the 
coefficient matrix si — A which is usually of large order. It also had conditioning 
problems which balancing ameliorated, but did not totally eliminate. 

In the new algorithm, the coefficient matrices involved in the solution of linear 
systems are A, which has the same order as the controller, and, for / = !.•••. p, 
sI- 2 —Ag, which is of order 2. Particularly when dealing with a large order structural 
model, the coefficient matrices used by the new algorithm are much smaller than 
the matrix si — A used by the old, so there is much less opportunity for round-off 
error. Any conditioning problems coming from the interaction of the frequency 
represented by s = juj and the z-th structural mode in the old method is isolated 
in the new method to calculating the denominator term u)'f — uj ' 2 + 2jQujjUj in 
equation (13); and this only gives numerical problems if uj is so close to jJ; that 
truncation occurs in forming the difference, and is so small that the (small) real 
part uj'f — uj ' 2 is a significant part of the whole term. 

Summary 

An efficient and novel procedure has been developed for the calculation of 


27 



the frequency response function of a large order, flexible system implemented 
with a linear, time invariant control system. The procedure takes advantage of 
the highly structured sparsity of the system matrices of the plant in normal mode 
coordinates. This reduces the computational cost from a quadratic function of the 
order of the system to a linear function, thereby permitting the practical frequency 
analysis of systems of much larger order than by traditional, full-matrix means. 
Formulations have been given for both open and closed loop systems. Numerical 
examples were presented wherein the advantages of the present formulation over 
traditional approaches, both in speed and in accuracy have been demonstrated. 
When exercised on the largest systems, the new Mex-code was about 40 times as 
fast as the old FORTRAN code when many frequency points were used while the 
advantage increased to a factor of about 80 or better when the calculation involved 
few frequency points. In this latter case, the new M-code was over 200 times as 
fast as the old M-code. 
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Table 1 Time (seconds) to Calculate Frequency Response Function 


System 

Size" 

Freq. 

vector 

length 

Feed- 

forward 

present 

New 

Mex-code 

New M-code 

Old 

FORTRAN 

code 

Old 

M-code 

1/1 

24+39 

200 

no 

1.98 

11.18 

1.26 

7.57 

1/1 

24+39 

200 

yes 

2.47 

13.50 

1.10 

7.62 

1/1 

24+39 

2000 

no 

24.78 

108.87 

10.08 

71.18 

1/1 

24+39 

2000 

yes 

19.60 

135.20 

10.36 

71.92 

3/5 

184+39 

200 

no 

5.03 

18.53 

29.84 

287.93 

3/5 

184+39 

200 

yes 

3.85 

14.18 

27.50 

290.25 

3/5 

184+39 

2000 

no 

46.08 

147.00 

213.69 

2818.78 

3/5 

184+39 

2000 

yes 

35.77 

188.93 

200.09 

2830.32 

10/27 

1406+39 

200 

no 

60.40 

231.25 

5714.94 

49846.62 

10/27 

1406+39 

200 

yes 

73.77 

227.40 

5756.99 

49199.20 

10/27 

1406+39 

2000 

no 

591.50 

1922.75 

24928.01 

- 

10/27 

1406+39 

2000 

yes 

615.73 

1897.22 

24929.14 

- 


a In "m/n" in the first line, "m" is the number of inputs and "n" is the number of outputs. 

In "m+n" in the second line, "m" is the number of structural states and "n" is the number of controller states. 
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